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INTRODUCTION 


Prediction  methods  of  multi-dimensional,  unsteady,  incompressible 
flows  can  now  investigate  many  engineering  problems  of  practical 
importance.  This  is  possible  principally  because  of  recent  advancements  in 
supercomputer  technology.  Besides  the  flow  being  unsteady,  the  problem 
geometry  is  often  complicated.  Airfoils,  intakes,  blade  passages,  heat 
exchangers,  recessed  cavities,  and  wall  mounted  bodies  are  just  a  few 
examples.  In  these  problems,  the  numerical  results  can  provide  sufficient 
detail  to  gain  a  comprehensive  understanding  of  the  basic  unsteady  flow 
physics. 

Basically,  three  distinct  approaches  are  available  for  solving  unsteady 
incompressible  flows.  One  particular  class  of  methods  is  founded  on  the 
fractional-step  procedure  (also  called  the  projection  method)  proposed  by 
Chorin  [1]  over  twenty-five  years  ago.  This  procedure  capitalizes  on 
Helmholtz's  decomposition  theorem  in  which  the  convection  and  diffusion 
components  of  the  unsteady  Navier-Stokes  equations  are  together 
decomposed  into  the  sum  of  a  divergence-free  vector  (velocity)  and  an 
irrotational  vector  (pressure  gradient).  Splitting  the  operators  in  this  way 
permits  the  introduction  of  an  intermediate  velocity  which  provides  the  only 
direct  link  between  the  pressure  and  physical  velocity  fields.  The  Navier- 
Stokes  equations  are  replaced  by  a  two-step  process  which  involves  solution 
of  an  unsteady  convection-diffusion  equation  for  the  intermediate  velocity 
followed  by  an  update  of  the  physical  velocity  field.  In  the  original  two- 
dimensional  scheme  by  Chorin,  the  method  of  artificial  compressibility  was 
implemented  to  establish  a  separate  equation  for  solution  of  the  pressure 
variable.  Given  the  intermediate  velocity  field,  the  physical  velocity  and 
pressure  fields  were  successively  updated  until  the  iterations  converged. 
Convergence  was  defined  by  the  allowable  error  in  the  incompressibility 
constraint.  At  each  time  increment,  many  iterations  were  necessary  to 
minimize  that  error.  Herein,  a  variation  of  Chorin’s  original  fractional -step 
method  is  utilized.  But  before  presenting  the  details  of  the  technique, 
alternate  methods  that  one  can  choose  for  solving  unsteady  incompressible 
flows  are  discussed. 

Full  predictor-corrector  type  schemes  have  also  been  applied  to 
unsteady  incompressible  flows.  One  approach  (referred  to  as  a  pressure- 
based  method)  involves  replacing  the  continuity  equation  with  a  pressure- 
Poisson  type  equation  which  is  derived  by  taking  the  divergence  of  the 
momentum  equation  and  invoking  the  incompressibility  constraint.  The 
pressure-Poisson  equation  along  with  momentum  are  converged  to  each  new 
time  level  by  a  semi-implicit  prediction  and  correction  sequence. 
Convergence  is  commonly  achieved  by  minimizing  the  residual  (or  error)  of 
the  discretized  pressure-Poisson  equation.  For  steady  flows,  large  time 
steps  are  allowed  to  accelerate  convergence  to  steady-state.  The  most 
recent  release  of  this  approach  is  called  PISO  (Pressure-Implicit  with 
Splitting  of  Operators)  [2].  Its  roots  lie  in  the  SIMPLE  (Semi-Implicit 
Method  for  Pressure-Linked  Equations)  algorithm  [3].  In  PISO,  the 
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predictor-corrector  steps  arise  from  a  splitting  of  operators  in  the 
discretized  momentum  and  pressure-Poisson  equations.  For  the  steady  flow 
problem.  SIMPLE  iterates  to  steady-state  convergence  while  PISO  advances 
in  time.  Thus,  PISO  is  equally  applicable  to  unsteady  and  transient  flows. 

Due  to  the  explicit  definitions  in  the  corrector  steps  however,  under¬ 
relaxation  may  be  necessary  in  PISO  to  sustain  convergence  during  the 
predictor-corrector  solution  sequence  [4], 

The  third  basic  approach  for  predicting  unsteady  incompressible  flows 
is  the  "method  of  pseudo-compressibility"  [51.  This  method  was  introduced 
by  Chorin  [6]  for  computing  steady  incompressible  flows.  Chorin's  idea 
resolved  the  difficulty  of  having  an  explicit  definition  for  the  pressure 
variable  in  the  incompressible  flow  equations.  He  proposed  to  couple  the 
continuity  and  momentum  equations  by  adding  a  pressure-like  time- 
dependent  term  to  the  continuity  equation.  This  term  physically  represents 
the  imposition  of  finite-speed  pressure  waves  onto  the  incompressible 
velocity  field.  The  pressure  waves  are  time-advanced  (along  with 
momentum)  until  they  theoretically  vanish.  At  convergence  the 
incompressibility  constraint  is  satisfied,  and  the  pressure  field  is  now  fully 
known.  For  the  steady  flow  problem,  it  is  important  to  note  that  the 
advancement  of  the  velocity  and  pressure  fields  to  steady-state  is 
accomplished  strictly  through  pseudo  time.  Consequently,  the  intermediate 
solutions  are  not  time  accurate.  Recently  however,  the  method  of  pseudo¬ 
compressibility  was  extended  to  unsteady  incompressible  flows  (see  Soh  and 
Goodrich  [7]  or  Rogers  and  Kwak  [8]).  The  new  system  of  equations  is 
essentially  the  same  as  the  original  set  except  for  a  second  time-dependent 
velocity  term  in  the  momentum  equation.  While  the  original  time- 
dependent  term  represents  real  time,  the  second  term  allows  for  marching 
in  pseudo  time.  Thus,  pseudo  steady-state  convergence  is  achieved  for  each 
physical  time  increment.  When  using  this  approach,  one  should  be 
concerned  about  the  local  error  tolerance  associated  with  pseudo  steady- 
state  convergence  relative  to  the  leading  temporal  truncation  error  of  the 
physical  system. 

Presented  here  is  a  solution  technique  for  solving  unsteady 
incompressible  flows  that  is  based  on  the  original  fractional-step  method 
developed  by  Chorin  along  with  the  variations  proposed  by  Kim  and  Moin  [9] 
and  Rai  and  Moin  [10].  The  present  technique  differs  from  those  previous 
works  in  three  basic  aspects.  First,  the  grid  is  semi-staggered,  meaning  that 
the  pressure  is  computed  at  the  grid  cell  centers  while  the  velocity 
components  are  collocated  with  the  grid  points.  Construction  of  the 
computational  molecule  in  this  fashion  greatly  facilitates  derivation  of  a 
consistent  set  of  boundary  conditions  for  the  intermediate  velocity.  Second, 
strong  coupling  is  maintained  between  the  pressure  and  velocity 
components  through  introduction  of  a  fourth-order-accurate  compact 
differencing  scheme  for  computing  the  pressure  gradients.  Lastly,  the 
modified  strongly  implicit  procedure  is  implemented  for  solution  of  the 
residual  form  of  the  discretized  pressure-Poisson  equation  in  three 
dimensions.  For  the  moment,  the  temporal  and  spatial  discretization  of  the 
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governing  equations,  implementation  of  the  solution  algorithms  and  test 
cases  are  presented  in  Cartesian  coordinates.  The  numerical  accuracy  of  the 
solution  technique  and  the  consistency  of  the  boundary  conditions  are 
verified  by  simulating  an  exact  solution  of  the  two-dimensional  Navier- 
Stokes  equations.  The  technique  is  then  used  to  predict  the  flow  in  a  three- 
dimensional  shear-driven  cavity  at  Reynolds  numbers  of  2000  and  3200.  At 
these  Reynolds  numbers  the  flow  is  unsteady  and  strictly  laminar.  The  test 
cases  include  comparisons  between  the  numerical  results  and  published 
experimental  data. 


FORMULATION  AND  DISCRETIZATION 

The  mathematical  system  which  governs  unsteady  incompressible  flows 
is  the  Navier-Stokes  (NS)  equations  with  continuity.  In  index  notation,  the 
conservative  non-dimensional  form  of  this  primitive-variable  system  in 
three-dimensional  (3D)  Cartesian  coordinates  is 


5uf  SUjU.  3p  l  32u, 

- -  H - -  = - —  -I - 

dt  dXj  dxi  Re  dxjdxj  ’ 


dUj 

dx, 


=  0, 


(la) 

(lb) 


where  the  velocity  vector  ui  =  <  u,v,w  >T  and  p  is  the  pressure.  The 
Reynolds  number  Re  =  UL/  v;  v  is  the  kinematic  viscosity  and  U  and  L  are 
the  reference  velocity  and  length,  respectively. 

Here,  Chorin's  original  fractional-step  method  along  with  the  variations 
proposed  by  Kim  and  Moin  [9]  and  Rai  and  Moin  [10]  are  combined  to 
temporally  discretize  the  system  in  equation  (1).  Kim  and  Moin  time-split 
the  convective  and  diffusive  terms  in  the  momentum  equation  using  an 
explicit  second-order-accurate  Adams-Bashforth  method  and  an  implicit 
Crank-Nicolson  scheme,  respectively.  Use  of  the  Crank-Nicolson  scheme 
eliminates  the  viscous  stability  restriction  which  can  be  severe  near  the  wall 
even  in  high-Re  flows.  As  pointed  out  by  Rai  and  Moin  however,  the  Adams- 
Bashforth  method  is  unconditionally  unstable  when  applied  to  the  linear 
convection  equation.  In  their  solution  technique,  they  chose  to  time-split 
the  convective  terms  using  a  third-order-accurate  Runge-Kutta  procedure. 
The  procedure  was  subsequently  improved  by  Le  and  Moin  [11]  in  terms  of  a 
reduced  computational  requirement  by  solving  the  pressure-Poisson 
equation  only  in  the  final  step  of  the  Runge-Kutta  procedure.  For  the  linear 
convection  equation,  this  particular  procedure  is  conditionally  stable. 
Contrasting  the  numerical  stability  of  the  Adams-Bashforth  method  and 
Runge-Kutta  procedure  is  presented  later  in  the  Numerical  Stability  section. 

The  Runge-Kutta/ Crank-Nicolson  solution  sequence  adopted  here  for 
time-advancement  of  the  system  of  equations  in  (1)  involves  three  complete 
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steps  to  update  the  physical  velocity  components  and  pressure  variable. 
This  sequence  has  the  form 
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The  parameter  n  denotes  the  particular  time  level  and  m  signifies  the 
three  steps  of  the  Runge-Kutta  procedure;  the  velocities  u”  -1  =  0,  u,n0  =  u” 

and  u,n  3  =  u,n+1 .  A  pressure  variable  0  replaces  the  actual  pressure  p  in  the 
velocity  update  equation  due  to  the  implicit  treatment  of  the  diffusion  term. 
An  exact  relationship  between  0  and  p  is  given  by  the  simple  expression  in 
(2f).  As  in  the  pressure-based  method,  a  Poisson  type  equation  is  derived 
for  solution  of  the  pressure  variable  0  by  taking  the  divergence  of  the 
velocity  update  equation  in  (2b)  and  enforcing  continuity  [equation  (2c)]. 
The  result  shown  in  equation  (2d)  provides  for  the  solution  of  0  in  terms  of 
the  divergence  of  the  intermediate  velocity  field  up  This  velocity  field  does 
not  satisfy  incompressibility  because  the  pressure  gradient  is  eliminated 
from  the  momentum  equation.  Since  the  intermediate  velocity  and 
pressure-Poisson  equations  are  solved  independently,  it  is  not  necessary  to 
treat  the  Runge-Kutta  coefficients  explicitly  for  updating  the  physical 
velocity  field.  A  second  scalar  field  ®  is  therefore  introduced  of  which  the 
pressure  can  be  computed  exactly  anytime  during  the  computation  through 
expressions  (2e)  and  (2f).  It  should  be  noted  that  by  retaining  the  Crank- 
Nicolson  scheme  for  the  diffusive  terms,  the  solution  accuracy  of  this 
fractional-step  technique  is  second  order  in  time. 
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One  should  notice  that  the  system  of  equations  in  (2)  is  now  in  a 
nonconservative  form.  This  form  facilitates  implementation  of  local  high- 
order  spatial  accuracy  for  the  convective  terms.  Inasmuch  as  flow 
discontinuities  are  absent  from  the  applications  considered  here,  the 
nonconservative  form  should  give  results  comparable  to  that  of  the 
conservative  form  [10].  Except  for  the  convective  terms,  all  spatial 
derivatives  shown  in  the  Runge-Kutta/ Crank  Nicolson  solution  sequence  are 
discretized  by  second-order  central  differences.  The  convective  terms  are 
explicit  and  are  spatially  discretized  by  a  third -order- accurate  five-point 
upwind-biased  stencil.  Velocity  points  which  reside  outside  the  grid 
boundaries  in  the  stencil  are  determined  by  extrapolating  the  field  results  to 
the  second  order  [12].  This  spatial  discretization  and  extrapolation  scheme 
leads  to  numerical  solutions  which  are  globally  second-order  accurate  in 
space. 


SOLUTION  METHODOLOGY 

In  this  section,  specific  details  are  presented  about  the  procedures  used 
for  solving  the  intermediate  velocity  equation  independent  of  the  pressure- 
Poisson  equation.  Establishing  proper  coupling  between  the  velocity 
components  and  pressure  variable  is  an  essential  ingredient  of  the  method 
to  inhibit  spurious  oscillations.  As  illustrated  in  detail  by  Patankar  [13],  fully 
staggered  grids  satisfy  this  demand  ideally,  but  require  extrapolation  or 
reflection  of  the  field  results  for  assigning  values  to  the  point  velocities 
which  lie  outside  the  geometric  boundaries.  The  extrapolation  becomes 
particularly  cumbersome  for  spatial  differencing  approximations  that  are 
higher  than  second  order.  One  can  also  implement  a  semi-staggered 
pattern  in  which  the  velocity  components  are  positioned  at  the  grid  points 
and  the  pressures  are  located  at  the  grid  cell  centers.  Besides  having  the 
velocity  field  now  defined  directly  on  the  boundary,  this  grid  also  offers  the 
advantage  of  deriving  a  consistent  set  of  boundaiy  conditions  for  the 
intermediate  velocity.  Without  proper  care  in  the  discretization  definitions, 
solution  difficulties  with  this  choice  will  arise  because  the  grid  cell 
pressures  can  actually  become  uncoupled  from  their  adjacent  neighbors  [14]. 
Described  below  is  a  computational  molecule  and  solution  procedure  that 
exploits  the  advantages  of  the  semi-staggered  grid  pattern  while  maintaining 
strong  coupling  between  the  velocity  components  and  pressure  variable. 

The  two-dimensional  (2D)  computational  molecule  in  Figure  1  shows 
the  relative  positions  of  the  velocity  components  and  the  pressure  variable 
on  the  grid.  The  velocity  components  are  collocated  with  the  grid  points, 
but  the  pressure  variable  is  staggered.  The  source  term  in  the  pressure- 
Poisson  equation  is  computed  by  averaging  the  grid  cell  corner  velocities.  If 
the  pressure  gradient  in  the  velocity  update  equation  is  now  determined  by 
averaging  the  pressures  at  the  grid  cell  centers,  then  the  semi -staggered 
grid  which  uncouples  the  cell  pressures  from  their  adjacent  neighbors  will 
be  recovered.  To  avert  this  dilemma,  the  pressure  gradients  are  first 
computed  at  the  cell  interfaces  (indicated  by  the  arrows  in  Figure  1)  by  a 
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NOTATION 
X:  VELOCITY  NODE 

(Uj,  Uj) 

O:  PRESSURE  CELL 
(P,  <b 

PRESSURE  GRADIENT 

(acv3x|) 


Figure  1.  Computational  Molecule  Used  for  Discretizing 
the  Incompressible  Navier-Stokes  Equations 


fourth-order-accurate  compact  differencing  formula.  The  pressure 
gradients  needed  to  update  the  point  velocities  are  obtained  by  averaging 
the  appropriate  values  at  the  cell  interfaces.  For  example,  the  u-velocity  at 
point  5  is  updated  by  averaging  the  pressure  gradients  at  cell  interfaces 
marked  N  and  S.  Likewise,  the  pressure  gradients  at  the  W  and  E  cell 
interfaces  are  averaged  together  to  update  the  corresponding  v-velocity.  By 
computing  the  pressure  gradients  in  this  manner,  the  velocity  components 
and  pressure  variable  remain  strongly  coupled.  Also,  no  superfluous  errors 
are  introduced  because  the  accuracy  of  the  pressure  gradients  are  within  the 
leading  truncation  error  of  the  overall  solution  technique. 

After  solving  the  pressure  Poissdn  equation  for  Om,  the  pressure 
gradient  is  computed  by  the  fourth-order-accurate  compact  scheme.  For 
example,  at  the  cell  interface  labeled  S  in  Figure  1, 


f  k!  +  22f  j  +  fj_j 


(3) 


where  fi  =  OOm/5x)j.  The  index  i  symbolizes  the  center  point  at  cell 
interface  S  and  i  +  1/2  denotes  the  cell  center  marked  (b).  Likewise,  at- 
cell  interface  W 


+  22gj  +  gj_j 


(4) 


where  gj  =  OOm/5y)j.  The  interface  center  point  is  denoted  by  index  j  and 
j  +  1/2  signifies  cell  center  (d).  Application  of  the  compact  scheme  to  all 
cell  interfaces  in  the  domain  produces  an  algebraic  set  of  equations  which 
can  be  solved  by  a  standard  tridiagonal  solver.  The  required  boundary 
conditions  are  determined  from  the  respective  velocity  update  equation 
(2b).  For  example,  application  of  equation  (3)  at  cell  interface  i  =  S  in  Figure 
1  requires  a  boundary  definition  for  fi- 1 .  This  definition  is  supplied  by  the 
velocity  update  equation  as 


At 


'30  Y" 

<  3X  Jj-i 


(5) 


where  u”:™  and  u™,  are  the  averaged  velocity  boundary  values  centered 
between  grid  points  1  and  4.  The  pressure  gradient  gj-i  in  equation  (4)  is 
treated  similarly.  There 


(vjf -v“  J 
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where  v":™  and  v™.j  are  the  averaged  velocity  boundary  values  centered 

between  grid  points  1  and  2.  Finally,  one  should  note  that  the  pressure 
gradients  for  updating  the  w-velocity  components  are  computed  in  an 
analogous  manner. 

The  intermediate  velocity  equation  is  solved  by  using  the  approximate- 
factorization  technique.  Because  this  equation  is  not  coupled  explicitly  to 
the  pressure-Poisson  equation  in  terms  of  the  scalar  Om,  it  is  solved  only 
once  for  each  intermediate  velocity  component  at  the  outset  of  each  Runge- 
Kutta  step.  Thus,  the  computational  efficiency  of  this  fractional-step 
technique  is  governed  primarily  by  the  CPU  time  required  to  solve  the 
pressure-Poisson  equation.  This  equation  satisfies  the  incompressible 
constraint  to  within  a  user-specified  error  tolerance,  and  its  rapid 
convergence  by  an  iterative  scheme  is  essential  to  the  overall  solution 
procedure.  For  this  purpose  the  modified  strongly  implicit  (MSI)  scheme, 
developed  by  Schneider  and  Zedan  [15],  is  implemented  in  a  residual  form 
of  the  pressure-Poisson  equation.  The  computational  efficiency  of  the  MSI 
scheme  over  point-successive  relaxation  and  alternating  directional  implicit 
techniques  was  demonstrated  by  Jordan  and  Spaulding  [16]  for  2D  grid- 
generation  and  by  Jordan  [17]  for  solving  the  vorticity-stream  function 
equations  in  2D  steady-state  flow  problems.  Extension  of  this  solver  to  the 
3D  pressure-Poisson  equation  is  a  straightforward  process.  It  can  be  viewed 
as  an  implicit  solver  of  a  planar  surface  (either  x-y,  x-z  or  y-z)  that  passes 
through  the  computational  volume  in  a  direction  (z,  y  or  x,  respectively) 
which  is  normal  to  the  surface.  The  2D  solver  sweeps  through  the  volume 
analogous  to  the  manner  an  implicit  line  solver  sweeps  along  a  planar 
surface. 

The  residual  form  of  the  discretized  pressure-Poisson  equation  appears 
as 


[A,  +A2  +Aa]5‘+1  =-R\  (7) 

where  the  coefficient  matrix  has  components  A(  =  52  /5x2,  the  increment 

-Om  <:  and  R'  is  the  residual.  The  parameter  t  denotes  successive 
updates  to  the  variable  Om  and  m  is  the  particular  Runge-Kiitta  step. 

For  the  pressure-Poisson  equation,  the  residual  is  defined  as 


d2®me  1  gu,m 
dxidxi  At  cxi 


(8) 


Solution  of  the  increment  5  in  equation  (7)  by  the  MSI  scheme  is  a  three- 
step  process  which  constitutes  a  single  iteration  t.  Each  step  reflects  a 
single  sweep  of  the  solution  scheme  through  the  computational  volume.  The 
directional  order  of  each  sweep  is  irrelevant  to  the  final  solution,  although 
the  number  of  total  iterations  required  for  convergence  may  sometimes  be 
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affected.  The  three-step  sequence  chosen  here  is  in  the  order  of  z,x,y  and 
has  the  form 

(A,  +  Aj)6^1  =  -Rf,  (9) 

where  i,j  are  cyclic  indices.  For  quickest  convergence,  the  residual  is 
recomputed  after  each  sweep.  The  full  3D  solution  is  now  dimensionally 
reduced  to  three  2D  implicit  solutions  for  which  the  MSI  scheme  is  directly 
applied.  Convergence  is  monitored  by  computing  the  root-mean-square 

(RMS)  of  the  residual.  At  convergence  R^g^O  meaning  that  the 
incompressibility  constraint  has  been  satisfied  to  within  the  error  tolerance 
specified  by  the  user. 

By  factoring  the  coefficient  matrix  in  this  way,  the  MSI  scheme  can  be 
implemented  easily.  This  scheme  was  originally  designed  by  Schneider  and 
Zedan  [15]  as  a  nine-point  implicit  solver  which  makes  it  well-suited  to 
handle  a  2D  curvilinear  form  of  the  pressure-Poisson  equation.  But  by 
conveniently  eliminating  the  appropriate  components,  it  can  be  applied 
effectively  to  the  2D  Cartesian  coordinate  system.  If  the  pressure-Poisson 
equation  is  recast  into  a  residual  form  first,  rather  than  implementing  the 
MSI  scheme  directly,  the  solution  methodology  gains  several  distinct 
advantages.  From  an  accuracy  standpoint,  the  residual  which  is  a  necessary 
computation  in  the  residual  form  is  a  justified  gauge  for  monitoring 
convergence.  Also,  from  a  programmer's  viewpoint,  implementation  of  the 
MSI  scheme  is  greatly  simplified.  This  is  because  the  residual  term  is  zero¬ 
valued  everywhere  on  the  grid  boundaries;  therefore  tracking  those 
boundaries  in  the  computation  is  not  necessary.  This  latter  advantage  is 
particularly  important  when  computing  flow  in  domains  having  internal 
boundaries. 


BOUNDARY  CONDITIONS 

The  impetus  for  collocating  the  velocity  components  with  the  grid 
points  is  to  facilitate  derivation  of  a  consistent  set  of  boundary  conditions  for 
the  intermediate  velocity  components  along  no-slip  walls.  For  the  semi- 
staggered  grid,  a  proper  definition  for  the  intermediate  velocity  along  no¬ 
slip  wall  boundaries  begins  by  projecting  equation  (2a)  onto  the  wall.  The 
result  is 


u  -  u 


n,m-l 


ymAt  a2(u,m  +u,n  m“I) 

Re  dxjdxj 


(10) 


which  is  implicit  in  u"1.  The  components  of  u"1  can  be  approximated  in 
terms  of  u(n~1  up  to  second -order  accuracy  in  time  by  noting  that 
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(11) 


un.m  _  =  _AtVOm"'  -O(At)2  =  -  0™“'  -O(At)2. 


Since  u"  m  =  un,m  1  along  the  wall,  equation  (11)  shows  that 

u™  =  u,m_1  +0(At)2.  By  substituting  this  result  into  equation  (10),  the 
definition  for  the  intermediate  velocity  becomes 
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-U 


n,m-l 


YmAt  a2(u;n-1  +  u"-m-1) 
Re 


+  O(At)3, 


(12) 


which  is  consistent  with  the  temporal  accuracy  of  the  field  solution.  To 
maintain  the  same  spatial  accuracy  along  the  boundary  as  in  the  field,  the 
diffusive  term  is  approximated  by  second-order-accurate  one-sided 
differences.  The  consistency  and  accuracy  of  this  boundary  condition  is 
tested  in  the  next  section. 


Two  important  concerns  are  satisfied  by  staggering  the  pressure 
variable  in  the  computational  molecule.  Besides  inhibiting  spurious 
oscillations  in  the  numerical  solutions  [13],  no  ad  hoc  wall  boundary 
conditions  are  necessary  for  the  pressure  variable.  This  can  be  illustrated 
effectively  by  the  following  simple  example.  The  residual  for  the  one¬ 
dimensional  form  of  the  discretized  pressure-Poisson  equation  at  the  center 
cell  point  i  =  a  in  Figure  1  is 


R'  =  («>,., -2<t,+<l>H1),',"-^[(us+u2)-(u4+u,)]".  (13) 

The  u-velocity  update  equation  applied  at  the  no-slip  wall  (u,n  =  0)  to  the  left 
of  grid  cell  (a)  becomes 


(i4) 

Substitution  of  this  equation  into  equation  (13)  in  terms  of  (u4  +u,)  reveals 
that  the  evaluation  of  the  residual  closest  to  the  wall  boundary  does  not 
require  an  ad  hoc  definition  for  the  pressure  variable  at  the  wall  (or  outside 
the  wall)  in  terms  of  the  field  values.  As  a  consequence,  the  pressure  field 
remains  conservative  and  second -order  accurate,  and  the  rate  of 
convergence  significantly  improves  because  an  ad  hoc  definition  commonly 
appears  as  a  Neuman  type  boundary  condition.  Extension  of  this  example  to 
three-dimensions  will  yield  the  same  result. 


SOLUTION  ACCURACY 

The  temporal  and  spatial  accuracy  of  the  overall  solution  procedure  as 
well  as  the  consistency  of  the  intermediate  velocity  boundary  conditions  can 
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be  verified  by  simulating  an  exact  solution  of  the  2D  unsteady  NS  equations. 
Through  reduction  of  the  mesh  spacing  under  a  constant  CFL  value,  the 
numerical  results  should  display  second-order  improvement  as  given  by  the 
leading  truncation  error.  The  exact  solution  has  the  form  [1] 


u(x,y,t)  =  -cos(x)sin(y)e~2t 

(15a) 

v(x,y,t)  =  sin(x)cos(y)e~2t 

(15b) 

p(x,y,t)  =  -)4[cos(2x)  +  cos(2y)]e"lt 

(15c) 

The  simulation  was  initialized  by  the  exact  solution  for  time  t  =  0.0.  The 
spatial  domain  was  defined  as  0  <  x,y  <  k.  Using  an  exact  set  of  boundary 
conditions  for  the  physical  velocity  with  the  definition  in  equation  (12)  for 
the  intermediate  velocity  (including  convection  where  appropriate),  the 
initial  flow  field  was  advanced  to  the  same  physical  time  for  each  grid  tested. 
In  this  case,  t  =  0.35  seconds,  which  corresponds  to  a  reduction  of  the 
maximum  boundary  velocity  to  one-half  of  its  initial  value.  The  results  are 
plotted  in  Figure  2  where  the  grid  factor  is  the  ratio  of  the  grid  spacing  in 
the  reference  grid  (11  x  1 1  uniform)  to  the  refined  grid.  An  exact  error  was 
determined  as  the  absolute  maximum  difference  between  the  numerical 
predictions  and  the  exact  solution  normalized  by  the  maximum  value  in  the 
domain.  Figure  2  clearly  shows  a  linear  reduction  of  the  velocity  and 
pressure  errors  by  refining  the  grid.  These  reductions  indeed  signify  a 
technique  which  is  second-order  accurate.  Thus,  the  numerical  accuracy  of 
the  solution  technique  is  verified  along  with  the  consistency  of  the  boundary 
conditions  for  the  intermediate  velocity. 


NUMERICAL  STABILITY 

Many  previous  applications  of  the  fractional-step  method  to  unsteady 
viscous  flows  temporally  discretized  the  convective  term  by  the  Adams- 
Bashforth  method  (see  Antonopoulos-Domis  [18]  and  Kim  and  Moin  [9]  for 
example).  This  scheme  holds  the  advantages  of  being  second -order- 
accurate,  one-step  and  low-storage  (only  information  from  the  previous  time 
step  needs  to  be  saved).  Unfortunately,  a  Fourier  stability  analysis  will  reveal 
that  the  method  is  unconditionally  unstable  if  applied  to  the  linear 
convection  equation.  As  shown  in  Figure  3  however,  this  instability  is  weak 
for  Courant  numbers  less  than  0.5  which  explains  why  the  method  has  been 
successful  for  solving  viscous  flow  problems.  But  for  certain  viscous  flow 
problems  where  at  some  point  Re— >-x  during  the  simulation,  for  example  in 
an  impulsively  started  flow,  very  small  CFL  values  would  be  necessary  (at  the 
outset  in  this  example)  to  sustain  stable  solutions. 

The  fractional-step  technique  developed  here  is  expected  to  posses 
strong  stability  characteristics.  As  pointed  out  earlier,  a  third -order- 
accurate  three-step  low-storage  Runge-Kutta  procedure  was  chosen  to  meet 
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Figure  2.  Improvement  in  the  Exact  Errors  of  Velocity  and  Pressure 
with  Grid  Refinement  for  a  Constant  Courant  Number 
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Figure  3.  Numerical  Stability  of  the  Adams-Bastiforth  Method 
Applied  to  the  Linear  Convection  Equation 
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this  need  in  lieu  of  the  Adams- Bashforth  method.  This  procedure  is 
conditionally  stable  when  applied  to  the  linear  convection  equation.  As 
shown  in  Figure  4,  when  the  procedure  is  combined  with  third-order- 
accurate  upwind-biased  spatial  differences,  the  linear  convection  equation  is 
always  stable  for  Courant  numbers  a  <  1.6.  In  the  case  of  the  impulsively 
started  flow  problem,  reasonable  time  increments  would  be  permitted  early 
in  the  simulation.  Obviously,  as  the  relative  importance  of  a  viscous  term 
increases,  the  stability  limit  shown  in  Figure  4  would  be  relaxed. 


RESULTS  AND  DISCUSSION 

In  this  section  we  present  the  simulation  results  of  a  3D  shear-driven 
cavity  flow  at  moderate  Reynolds  numbers  of  2000  and  3200.  At  these 
Reynolds  numbers  the  flow  is  unsteady  and  fully  laminar.  In  the  discussion 
of  the  numerical  results,  we  include  qualitative  and  quantitative  comparisons 
to  the  experimental  data  previously  published.  Rhee  et  al  [19]  published 
flow  visualization  results  of  the  cavity  flow  at  Re  =  2000  and  Re  =  3300  while 
Koseff  and  Street  [20,21]  presented  extensive  visualization  and  statistical 
flow  data  at  Re  =  3200.  In  those  experiments,  the  cavity  width  (x-direction) 
equaled  the  height  (y-direction).  The  length  of  the  span  (z-direction)  was 
varied.  The  visualization  evidence  showed  that  the  2D  recirculation  flow  (x-y 
planes)  can  be  characterized  basically  by  a  primary  vortex,  a  downstream 
secondary  eddy,  an  upstream  secondary  eddy  and  an  upper  secondary  eddy 
(see  Figure  5).  These  basic  features  have  been  predicted  many  times  by 
solution  techniques  that  solved  the  classic  2D  problem  (see  Ghia  et  al  [22] 
for  example).  In  the  spanwise  direction,  the  flow  forms  Taylor-Gortler-like 
(TGL)  vortex  pairs  and  a  lower  comer  vortex  at  the  spanwise  end -walls. 
These  characteristics  are  illustrated  in  Figure  6.  According  to  Koseff  and 
Street  [23],  the  impetus  manifesting  the  TGL  vortices  is  the  instability  of  the 
concave  free  shear  layer  that  separates  the  primaiy  vortex  from  the 
downstream  secondary  eddy.  These  vortices  form  just  above  the  concave 
surface  much  like  the  flow  experiments  between  rotating  cylinders  by  Taylor 
[24]  and  the  concave  boundary  layer  investigation  by  Gortler  [25].  At  the 
Reynolds  numbers  under  consideration,  the  visualization  results  showed  the 
TGL  vortices  meandering  slowly  along  the  cavity  bottom.  The  size  and 
number  of  vortices  depends  strongly  on  the  Reynolds  number.  However, 
Rhee  et  al  [19]  and  Koseff  and  Street  [20,22]  noted  that  the  spanwise  flow 
was  symmetric  about  the.  mid-span  plane  when  the  spanwise  aspect  ratio 
(SAR)  of  the  experimental  apparatus  was  extended  to  three.  Manifestation 
of  the  comer  vortex  is  a  consequence  of  the  shear  and  pressure  force 
adjustment  in  the  recirculating  flow  caused  by  the  no-slip  condition  along 
the  spanwise  end-wall.  Although  this  vortex  is  quasi-steady  at  these 
Reynolds  numbers,  its  size  and  extent  strongly  influence  the  TGL  vortices. 
Thus,  sufficient  grid  resolution  must  be  provided  for  both  the  TGL  vortices 
and  the  comer  vortex  to  properly  capture  the  characteristics  of  this  truly 
unsteady  flow. 
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Figure  4.  Numerical  Stability  of  the  Linear  Convection  Equation  Using  the  Third-Order-Accurate 
Three-Step  Runge-Kutta  Procedure  and  Third-Order-Accurate  Upwind-Biased  Spatial  Differences 


U  =  1  V  =  0 


Figure  5.  Sketch  of  the  Basic  Features  of  Recirculation  in  the 
Two-Dimensional  Shear-Driven  Cavity  Flow  Problem 


Figure  6.  Sketch  of  the  Taylor-Gortler-Like  Vortex  Pairs 
as  Observed  in  the  Flow  Visualization  Experiments 
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In  both  simulations  of  the  3D  shear-driven  cavity  flow  ,  the  geometry 
was  modeled  with  unit  width  and  unit  height  (see  Figure  7).  The  span 
length  was  1  and  1/2  times  that  of  the  width  (or  height)  with  one  boundary 
modeled  as  a  plane  of  symmetry  at  the  mid-span  and  the  other  a  solid  end- 
wall;  thus  as  in  the  experiments,  SAR  =  3.0.  The  cavity  lid  moved 
horizontally  at  unit  velocity  and  the  no-slip  condition  was  applied  to  all 
boundaries  except  the  mid -span  plane  which  was  treated  numerically  as  a 
plane  of  symmetry.  Each  simulation  was  initialized  by  an  impulsively  started 
lid.  Extremely  small  time  steps  were  not  necessary  early  in  the  computation 
which  was  due  to  the  use  of  the  Runge-Kutta  procedure.  The  results  labeled 
t  =  0.0  depict  solutions  that  had  been  time-advanced  until  the  transient 
effects  of  the  impulsively  started  lid  became  negligible.  After  completion  of 
the  early  time  steps,  only  a  few  iterations  of  the  pressure-Poisson  equation 
were  necessary  to  reach  incompressibility  at  each  new  time  level. 

For  the  cavity  flow  test  case  at  Re  =  2000,  a  41  x  41  x  51  uniform  grid 
(x,  y.  z  directions)  was  used  with  CFL  =  2.4.  A  high  resolution  in  the 
spanwise  direction  was  chosen  to  insure  prediction  of  the  complete  TGL 
vortex  structure.  Three  sets  of  snapshots  are  shown  in  Figure  8  of  the 
recirculation  and  spanwise  flow  at  non-dimensional  times  of 
t  =  0.0,  8.0  and  16.0.  The  recirculation  velocity  vectors  typify  the  flow  at  the 
mid-span  plane.  For  clarity,  they  are  normalized  with  respect  to  their  own 
magnitudes.  The  spanwise  velocity  vectors  show  the  flow  at  a  plane  where  x 
=  0.77.  The  relative  time  t  =  0.0  is  an  actual  time  of  24  seconds  after 
impulsive  start  of  the  cavity  lid.  At  this  initial  time,  three  small  TGL  vortex 
pairs  span  the  cavity  bottom  with  a  large  corner  vortex  located  at  the  lower 
end-wall.  The  corresponding  recirculation  flow  displays  strong  2D  features 
which  is  a  consequence  of  the  weak  spanwise  flow.  After  an  additional  eight 
time  units  however,  the  TGL  vortex  closest  to  the  mid-plane  nearly  doubled 
its  physical  size.  This  now  larger  vortex  pair  impacts  the  recirculation  flow 
locally  by  extracting  kinetic  energy  from  the  downstream  region; 
consequently  reducing  the  size  of  the  downstream  secondary  eddy.  At  time 
t  =  16.0,  three  TGL  vortices  clearly  appear  with  a  fourth  much  weaker  one 
positioned  near  the  corner  vortex.  The  prediction  of  four  TGL  vortices  at  Re 
=  2000  agrees  with  the  flow  visualization  results  reported  by  Rhee  el  al  [19]. 

In  the  flow  visualization  experiments  at  this  Reynolds  number,  the  TGL 
vortices  meandered  slowly  along  the  cavity  bottom  as  well  as  varied  their 
physical  size.  This  behavior  is  also  illustrated  in  Figure  9  where  five 
snapshots  of  the  computations  are  shown  at  unit  time  intervals  beginning 
with  t  =  26.  The  sequence  of  snapshots  show  creation  and  stationary  growth 
of  a  TGL  vortex  pair  directly  next  to  the  mid-span  plane.  Conversely,  the 
TGL  vortex  closest  to  the  end-wall  is  reduced  and  reverses  its  lateral  path 
three  times.  Furthermore,  this  particular  vortex  meanders  across 
approximately  15  percent  of  the  cavity  bottom.  The  remaining  two  TGL 
vortices  also  display  large  variations  in  their  physical  size,  but  meander 
comparatively  to  a  much  lesser  degree. 
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t  =  28 


t  =  29 


Figure  9.  Snapshots  of  Spanwise  Flow  Showing  Meander  and 
Physical  Size  ofTGL  Vortices  for  Re  =  2000  Test  Case 

20 


To  graphically  verify  the  claim  of  strong  coupling  between  the  velocity 
components  and  pressure  variable,  three  pressure  variable  profiles  are 
displayed  in  Figure  10.  The  pressure  variable  is  plotted  instead  of  the  actual 
pressure  to  avoid  the  natural  smoothing  of  the  diffusive  term  as  given  in 
equation  (2d).  The  profiles  represent  cuts  at  y  =  0.1  through  the  free  shear 
layers  that  separate  the  primary  vortex  from  the  secondary  eddies.  These 
profiles  extend  from  the  upstream  wall  (x  =  0.0)  to  the  downstream  wall  (x 
=  1.0)  at  the  mid-,  l/3rd  and  2/3rds  planes  along  the  cavity  span.  Notice 
that  the  profiles  depict  smooth  continuous  curves  that  are  devoid  of  spurious 
oscillations.  These  characteristics  illustrate  proper  coupling  between  the 
pressure  variable  and  velocity  components  as  given  by  the  discretization 
molecule. 

A51  x51  x  65  uniform  grid  was  selected  for  the  test  case  at  Re  =  3200 
and  a  CFL  value  of  1.5.  This  mesh  resolution  was  based  on  the  simulation 
reported  by  Freitas  et  al  [12]  where  a  32  x  32  x  45  grid  was  used  that  was 
non-uniform  in  the  recirculation  planes.  A  lower  CFL  value  was  chosen  for 
this  test  case  to  insure  capture  of  the  unsteady  flow  physics  in  the  cavity. 
Figures  11a  and  lib  show  typical  snapshots  of  the  unsteady  flow  (t  =  12.0). 
The  recirculation  velocity  vectors  (Figure  11a)  represent  the  flow  at  the 

mid-span  plane  whereas  the  spanwise  flow  (Figure  1  lb)  is  shown  at  the  x  = 
0.77  plane.  Once  again  for  clarity,  the  recirculation  velocity  vectors  are 
normalized  with  respect  to  their  own  magnitudes.  Four  TGL  vortex  pairs 
span  the  cavity  bottom.  Their  position  is  marked  by  grid  lines  referenced  to 
the  cavity  mid -span.  This  number  of  TGL  vortices  agrees  with  the 
experimental  observations  reported  by  Rhee  et  al  [19].  The  streamwise 
extent  of  each  vortex  pair  in  the  form  of  x-vorticity  contours  is  plotted  in 
Figure  12.  The  marking  of  grid  lines  in  that  figure  and  in  the  spanwise 
velocity  vectors  (Figure  lib)  help  quantify  the  intensity  of  each  vortex  with 
respect  to  its  spanwise  size  and  streamwise  extent.  Generally,  the 
streamwise  stretch  of  each  vortex  is  directly  dependent  on  its  size  and 
strength.  As  time  passes  however,  new  TGL  vortices  are  generated  which 
gain  physical  strength  and  size.  Also,  the  corner  vortex  over  the  same  time 
period  behaves  in  a  similar  manner.  For  example,  the  snapshots  in  Figure 
13  (t=75.0)  shows  a  large  TGL  vortex  straddling  the  mid-span  with  three 
smaller  ones  along  the  spanwise  plane.  The  larger  TGL  vortex  clearly 
dominates  the  local  downstream  region  of  the  cavity  near  the  mid-span  by 
completely  eliminating  development  of  the  secondary  eddy.  Normalized 
recirculation  velocity  vectors  representing  18  minute  sample  averages  at  the 
l/3rd  and  2/3rds  planes  from  the  spanwise  end-wall  are  shown  in  Figures 
13a  and  13b„  respectively.  At  these  planes,  downstream  and  upstream 
secondary  eddies  are  clearly  visible,  but  not  an  upper  secondary  eddy. 

Freitas  et  al  [12]  also  reported  this  result  and  both  results  agree  with  the 
flow  visualization  data  [19,20],  Furthermore,  the  l/3rd  plane  vectors  show 
the  primary  vortex  core  positioned  in  the  upper  right  quadrant  of  the 
recirculation  plane  whereas  the  core  at  the  2/3rds  plane  lies  close  to  the 
geometric  center.  This  result  also  agrees  qualitatively  with  the 
experimental  observations  and  is  due  to  the  proximity  influence  of  the  end- 
wall  on  the  recirculation  flow. 
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PRESSURE  VARIABLE 


Figure  10.  Pressure  Variable  Profiles  at  Y  =  0.1  from  the  Upstream  (x  =  0.0)  to  the  Downstream 
Wall  (x  =  1.0)  at  the  Mid ,  l/3rd  and  2/3 rds  Planes  in  the  Z  Direction 
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(a) 


(b) 


Figure  11.  Snapshots  of  the  Recirculation  (a)  and  Spanwise  Flow  (b)  in  the  Three-Dimensional  Shear- 
Driven  Cavity  at  Time  t=12.0  for  Re  =  3200  and  a  51  x  51  x  65  Uniform  Grid.  Recirculation  Vectors 
are  Shown  at  the  Cavity  Mid-Span  Whereas  the  Spanwise  Flow  is  at  Plane  X  =  0.77 


(b) 


Figure  13.  Normalized  Recirculation  Velocity  Vectors  Representing 
18-Minute  Sample  Averages  of  Shear-Driven  Cavity  Flow  at  Re  =  3200; 
(a)  l/3rd  and  (b)  2l3rds  Planes  from  the  Cavity  Spanwise  End-Wall 
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The  experimental  and  computed  time-averaged  u  and  v  velocity  profiles 
through  the  cavity  center  at  the  mid-span  are  compared  in  Figures  14(a)  and 
14(b),  respectively.  The  experimental  profiles  [20]  represent  five  minutes 
of  time-averaged  measurements,  whereas  the  present  computations  are 
averaged  over  18  minutes.  Also  included  in  the  figure  are  the  steady  2D 
velocity  profiles  for  a  51  x  51  uniform  grid  and  the  three  minute  time- 
averaged  results  of  Freitas  et  al  [12],  Comparisons  between  the  3D 
computations  and  experimental  results  are  quite  good.  On  the  other  hand, 
comparisons  between  the  2D  computation  and  both  3D  results  are  quite 
poor.  This  confirms  the  conclusion  drawn  by  Koseff  and  Street  [2Q],  that 
the  three-dimensional  effects  on  the  recirculation  flow  manifest  significant 
differences  between  the  2D  and  3D  mean  velocity  profiles.  One  should  note 
that  the  present  time -averaged  velocities  away  from  the  cavity  walls  agree 
better  with  the  experimental  data  than  the  corresponding  results  of  Freitas 
et  al  [12]  for  two  possible  reasons.  The  first  is  their  lower  sampling  window, 
but  more  importantly,  they  sacrificed  the  grid  field  resolution  by  clustering 
lines  near  the  cavity  walls  in  the  recirculation  planes. 

As  a  final  note,  after  completion  of  the  transient  phase  of  the 
computation  for  this  test  case,  the  CPU  requirement  was  approximately 

1.3xlOr5s  per  grid  point  per  Runge-Kutta  step  on  a  CRAY-YMP  platform. 

This  requirement  compares  competitively  to  the  computational 
requirements  reported  by  Kim  and  Benson  [26];  1.7xl0~5s  (iterative  time 
marching),  2.5xl0"5s  (simplified  marker-and-cell),  and  4.5xl0~5s  (PISO). 


CONCLUSIONS 

Details  of  a  new  fractional-step  technique  developed  specifically  for 
predicting  three-dimensional  unsteady  incompressible  flows  was  presented. 
Use  of  a  semi-staggered  grid  pattern  permitted  easy  derivation  of  a 
consistent  set  of  boundary  conditions  for  the  intermediate  velocity 
components  without  necessitating  interpolation  or  extrapolation  of  the  field 
data.  Also,  elimination  of  an  ad  hoc  boundary  condition  for  the  pressure 
solution  was  demonstrated.  The  difficulty  of  insuring  proper  local  coupling 
between  the  pressure  and  velocity  components  in  the  computational 
molecule  was  satisfied  by  computing  the  pressure  gradient  in  the  velocity 
update  equation  at  the  grid  cell  interfaces  using  fourth-order-accurate 
compact  differences.  Extension  of  the  modified  strongly  implicit  scheme 
for  solution  of  the  three-dimensional  pressure-Poisson  equation  in  residual 
form  provided  an  efficient  iterative  means  of  quickly  achieving 
incompressibility  at  each  sub-step  of  the  Runge-Kutta  procedure.  Both  the 
consistency  of  the  boundary  conditions  and  global  solution  accuracy  (second- 
order  in  time  and  space)  were  verified.  The  technique  gave  accurate 
qualitative  and  quantitative  results  of  the  three-dimensional  shear-driven 
cavity  flow  at  moderate  Reynolds  numbers  of  2000  and  3200.  The 
technique  is  also  computationally  competitive  against  other  schemes  of  the 
same  class. 
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CAVITY  WIDTH 


(b) 

Figure  1 4.  Comparison  Between  Experimental  and  Computational  Time- 
Averaged  Velocity  Profiles  Through  Cavity  Center  at  Mid-Span  for 
Re=3200.  Experimental  Data  Taken  from  Koseff  and  Street  [21]  and 
Reference  Simulation  from  Freitas  et  al  [12]: 

(a)  u  -  velocity,  (b)  v  -  velocity 
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